Chaotic dynamics around astrophysical objects with nonisotropic stresses 
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The existence of chaotic behavior for the geodesies of the test particles orbiting compact objects 
is a subject of much current research. Some years ago, Gueron and Letelier [Phys. Rev. E 66, 
046611 (2002)] reported the existence of chaotic behavior for the geodesies of the test particles 
orbiting compact objects like black holes induced by specific values of the quadrupolar deformation 
of the source using as models the Erez-Rosen solution and the Kerr black hole deformed by an 
internal multipole term. In this work, we are interesting in the study of the dynamic behavior 
of geodesies around astrophysical objects with intrinsic quadrupolar deformation or nonisotropic 
stresses, which induces nonvanishing quadrupolar deformation for the nonrotating limit. For our 
purpose, we use the Tomimatsu-Sato spacetime [Phys. Rev. Lett. 29 1344 (1972)] and its arbitrary 
deformed generalization obtained as the particular vacuum case of the five parametric solution of 
Manko et al [Phys. Rev. D 62, 044048 (2000)], characterizing the geodesic dynamics throughout the 
Poincare sections method. In contrast to the results by Gueron and Letelier we find chaotic motion 
for oblate deformations instead of prolate deformations. It opens the possibility that the particles 
forming the accretion disk around a large variety of different astrophysical bodies (nonprolate, e.g., 
neutron stars) could exhibit chaotic dynamics. We also conjecture that the existence of an arbitrary 
deformation parameter is necessary for the existence of chaotic dynamics. 

PACS numbers: 95.10.Fh, 05.45-a, 04.20Jb 



I. INTRODUCTION 



Although it is very usual in the literature that the stress 
tensor of several astrophysics objects matter -neutron 
stars, exotic stars- is approximated by that of a per- 
fect fluid, viz. nonviscous medium of total energy den- 
sity £ (mass density p = £/c 2 ), in which all stresses are 
zero except for an isotropic pressure P. We have to no- 
tice that in astrophysics, the most common object are 
that which present nonvanishing quadrupolar deforma- 
tion for the nonrotating limit, i.e., objects not deformed 
by current mass -rotation- either by arbitrary multipoles 
-arbitrary mass quadrupole or mass octupole- but intrin- 
sically deformed. Besides, it is clear that the approxima- 
tions is based on the fact that the shear stresses, e.g., 
those produced by elastic strain in the solid crust or by 
strong magnetic field, are generally negligible compared 
to the pressure, but it is also clear that when the pres- 
sure is assumed isotropic we refuse the possibility of in- 
trinsic deformation because it could obey to nonisotropic 
stresses during its formation process. 
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With the aim to follow the usual premise that the ex- 
terior gravitational field of the relevant astrophysical ob- 
jects can be modelled by stationary axially symmetric ex- 
act solutions to the Einstein-Maxwell field equations (see 
e -g- [2 for the case of neutron stars) we adopt space- 
times endowed naturally in the Weyl-Lewis-Papapetrou 
line element Q. This metric admits two Killing vectors, 
one timelike and other spacelike, which correspond to the 
desiderated symmetries mentioned above. In the litera- 
ture, we can find several exact solutions to the stationary 
axis-symmetric Einstein-Maxwell system equations (see 
[H for some of them), but just a few of them described 
objects intrinsically deformed. The most known are the 
corresponding members to the Tomimatsu-Sato family 
[f| which quadrupole deformation in the nonrotating 
limit is Q = —^-m 3 , being m — ^j 2 - where Mq is the 
mass monopole of the source and S a dimensionless pa- 
rameter, taking the values 6 = 1 for the Kerr solution 
and S = 2 for the Tomimatsu-Sato 5 = 2 solution. 

In order to clarify the physical interpretation of the 
Tomimatsu-Sato family (henceforth TS-family), we 
present briefly in the next paragraphs some interpreta- 
tions given to it. 

In the early 1970s, when the solution was obtained, 
Tomimatsu & Sato 0] and after Tanabe [1] showed that 
TS-family represents the gravitational field of rotating 
masses with angular momentum about the z-axis. But, 
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soon was realized that TS-family suffers of naked cur- 
vatures singularity. Papadopoulos and Xanthopoulos Q 
tried to resurrect interest in these solutions putting them 
in other context, they modified slightly the TS-metrics by 
a suitable analytic continuation and concluded that these 
solutions represent cylindrical symmetric spacetime, in- 
terpreting this time the solution as a beamlike-shaped 
pulse of gravitational radiation scattered by a cosmic 
string. Kodama & Hikida [l(| showed that the two points 
in the Weyl coordinates, which have been recognized as 
the directional singularities, are really two-dimensional 
surfaces and that these surfaces are horizons. They also 
showed that each of the two horizons has the topology of 
a sphere and concluded that this may indicates that TS- 
family describes the spacetime surrounding a new possi- 
bility of final states of gravitational collapse. 

In 2002, using the Sibgatullin's integral method [ill ]. 
Mielke, Manko and Sanabria-Gomez generalized a mem- 
ber of the vacuum TS-family, the Tomimatsu-Sato 6 = 2 
solution (henceforth TS2), to a most general spacetime 
in the electrovacuum case [6| possessing parameters for 
the electric charge Q, magnetic dipole p, and arbitrary 
quadrupolar deformation b. Besides, Berti & Stergioulas 
[l2| and Berti, White, Maniopoulou & Bruni studied 
this solution as possible model to describe the gravita- 
tional field of a Rapidly Rotating Neutron Stars (hence- 
forth RRNS). In that work the generalized version of the 
TS2 solution [6| was considered and it was concluded that 
this solution is a suitable model for the field of a neutron 
star but only in the faster rotation regime. The reason 
why the solution by Manko et al. cannot be used to 
model neutron stars close to the nonrotating limit, (i.e., 
when the source is slightly deformed -quasispherical) can 
be understood if we take into account that in this limit, 
neglecting the high order multipoles, [f| reduces to TS2 
leaving a nonvanishing deformation given by Q = —jm 3 , 
which differs substantially of the spherical symmetry due 
to the high mass of the star. 

In this point, following the concluding remarks by Ko- 
dama et al. fioj and Berti et al. (T3 |. we assume that 
TS-family can be used to describe the topology of the 
spacetime around relevant astrophysical objects like neu- 
tron stars, strange quarks or any other exotic final states 
of the gravitational collapse. This last statement and 
the fact that the members of the TS-family have non- 
vanishing quadupolar deformation constitute the initial 
point for our work because it enables us to consider the 
members of the TS-family as analytic closed form solu- 
tions for the gravitational field of relevant astrophysical 
objects with anisotropic stresses tensor. 

After having an analytic closed form solution for the 
gravitational field of a source, one of the interesting topic 
to study is the motion of the particles orbiting this source. 
The construction of astrophysical models which are able 



to give us a complete description of realistic astronomical 
systems has to take into account the behavior of the sur- 
rounding matter in order to compare with astronomical 
observation and emit a proposition about the validity of 
the model. In this work, we are interested in the nature 
of the dynamics -chaotic or regular- of the test particles 
which orbit two specific members of the TS-family, the 
TS2 and the solution by Manko et al. in the vacuum 
case, viz. Q = p, = 0. 

In general relativity the study of stochastic motions in 
deterministic system -deterministic chaos- has followed 
two main branches. The first one is the study of the 
geodesic motion of test particles in a given gravitational 
field (Bombclli and Calzetta Vieira and Letelier 

[IHO, Gueron and Letelier [l7l] and references therein). 
The other branch is the time evolution of the gravita- 
tional field itself (Motter and Letelier [llj], Hobill, Burd 
and Coley |l9(), which is relevant in cosmology (see e.g. 
Motter |2(|). This work is in line with the first scheme 
and is organized as follows. In Sec. II the TS2 solution 
and the Manko et aV solution are presented and a brief 
discussion of their features is given. In Sec. Ill the dy- 
namics of geodesies of test particles is analyzed. Finally, 
in Sec. IV a brief discussion about the obtained results 
is presented. 



II. THE PARTICULAR SPACE-TIMES 



A. Case I: The Tomimatsu— Sato 5 = 2 solution 



The most simple form of the metric for a stationary ax- 
isymmetric space time was given by Papapetrou [3[ and 
it can be written as 



ds 2 = f(dt — ivdip) 2 — f 



(dp 2 + dz 2 )+p 2 dcp% (1) 



here /, u> and 7 are functions of the quasicylindrical 
Weyl-Lewis-Papapetrou coordinates (p,z). The Weyl- 
Lewis-Papapetrou coordinates are related with the pro- 
late spheroidal coordinates (x, y) by means of the trans- 
formation 



P 



k 2 {x 2 -l){l-y 2 ). 



kxy. 



(2) 



with x > 1 and — 1 < y < 1 , then the metric is rewritten 
as 



ds 2 = f(dt- u<kpf - k 2 /- 1 e 2 ~< (x 2 - y 2 ) ( 



dx 2 



+ 



dy 2 



1 



+ (x 2 - 1)(1 - y 2 )d<p 2 



(3) 



The TS metrics were obtained from a solution to the 
Ernst equation in the vacuum case [2l| . which is given 
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by 

p 2 x 4 + q 2 y 4 — 2ipqxy(x 2 — y 2 ) — 1 . . 

^ = 2px(x 2 - 1) - 2iqy(l - y 2 ) ' ^' 

where p = (1 — q 2 ) 1 / 2 , q = J/m 2 and k = rap/5. The 5 
parameter is dimensionless, taking the values 5 = 1 and 
6 = 2 for the Kerr and the TS2 solutions, respectively. 
The metric functions derived from are [5| 

A _ 2mgC{\-y 2 ) 2l _ A 
1 B ' U A ' 6 p 4 (x 2 - y 2 Y ' 

(5) 

where 

A = p\x 2 -l) 4 + q\l-y 2 ) 4 -2p 2 q 2 (x 2 -l)(l-y 2 ) 

x {2{x 2 - I) 2 + 2(1 - y 2 ) 2 + 3(x 2 - 1)(1 - y 2 )} , 
S = {p 2 (x 2 + l)(x 2 -l)- <Z 2 (y 2 + l)(l-y 2 ) 

+ 2px(x 2 - l)} 2 + Aq 2 y 2 {px(x 2 - 1) 

+ (pa; + l)(l-y 2 )} 2 , 
C = -p 3 x(x 2 -l){2(x 2 + l)(x 2 ~l) + (x 2 + 3)(l-y 2 )} 

- p 2 (x 2 - l){4x 2 (x 2 - 1) + (3x 2 + 1)(1 - y 2 )} 

+ q 2 {px + - y 2 ) 3 . 

The physical sense of the parameters m and J is derived 
form the Simon's multipole moments [22| . The first three 
of the relativistic multipole moments calculated from (j4]) 
with the aid of the Hoenselaers-Perjes procedure [23| are 

The terms Pq and P2 denote the monopole and the 
quadrupole of mass, respectively. They are related to 
the mass and the deformation of the source. On the 
other hand, the term P± describes the angular momen- 
tum. From ([6]) is obvious that Pq, P2, and Pi are deter- 
mined only by two parameters, m and J, and means that 
the source only has mass and arbitrary angular momen- 
tum. 

Defining the parameter j as J/m 2 the quadrupole of 
mass for the case 6 = 2 can be written as P2 = Q = 
— 0.25to 3 (1 + 3j 2 ) . From the previous expression it is 
clear that the quadrupole deformation of TS2 solution al- 
ways is negative. For such reason we can affirm that the 
TS2 solution describes the spacetime around an oblate 
source. 



momentum (a = J/m), Q the electric charge, b a param- 
eter related with the mass quadrupole moment, and fi a 
parameter related with the dipolar magnetic moment. In 
addition, it has a remarkable feature, which is that the 
quadrupolar deformation, 

Q = m{5 -d-a{a-b)), (7) 

with 

fi 2 — m 2 b 2 
m 2 -(a- b) 2 - Q 2 ' 

d = I[ m 2-(a-6) 2 -Q 2 ], (8) 

depends directly on electromagnetic parameter Q and /z. 
This means that a test particle "sees" a source with a dif- 
ferent quadrupolar deformation than the real one. This 
feature could implies that the electromagnetic field can 
induce chaos in the dynamic geodesic for uncharged par- 
ticles orbiting this source [24J. In this work, we are inter- 
ested in the influence of the real deformation of the mass 
and not in the effective deformation. For that reason, we 
will choose Q = and /1 = 0, then the metric functions 
take the form 

E 27 E = -(l~y 2 )F 

1 D' 6 16K 8 (a; 2 -y 2 ) 4 ' ^ E 

(9) 

with 

E = {4[k 2 (x 2 - 1) + 6(1 - y 2 )] 2 + (a - b) 
[(a - b)(d - 6) - m 2 b](l - y 2 ) 2 } 2 
-16«V - 1)(1 - y 2 ){(a - b)[n 2 (x 2 - y 2 ) 
+26y 2 }+m 2 by 2 } 2 , 

D = {4(k 2 x 2 - 6y 2 ) 2 + 2Kmx[2K 2 (x 2 - 1) 

+ (26 + ab- b 2 )(l - y 2 )] + (a - b)[(a -b)[d- 6) 
-m 2 b](y 4 - 1) - 4d 2 } 2 + A V 2 {2k 2 (x 2 - 1) 
[nx(a — b) — mb] — 2mb6(l — y 2 ) + [(a — b) 
(d-S)- m 2 b](2Kx + m)(l - y 2 )} 2 , 

F = 8n 2 (x 2 -l){(a-b)[K 2 (x 2 -y 2 ) + 25y 2 } 
+m 2 by 2 }{Kmx[(2Kx + m) 2 — 2y 2 (25 + ab 
-b 2 )-a 2 + b 2 ]-2y 2 (Ud 
-m 2 b 2 )} + {A[n 2 (x 2 - 1) + 8(1 - y 2 )] 2 
+ (a - b)[(a - b)(d -5)- m 2 b](l - y 2 ) 2 } 
(1(2nmbx + 2m 2 6)[K 2 2 - 1) + 5(1 - y 2 )] 
+ (1 - y 2 ){(a - b)(m 2 b 2 - A5d) - (Anmx + 2m 2 
)[(a-b)(d-5)-m 2 b}}). (10) 



B. Case II: Solution by Manko, Mielke and 
Sanabria— Gomez 



This solution has five relevant and independent param- 
eters: m the gravitational mass, a the specific angular 



j^(r+ — r_) with r± = 
If, in addition, we put 



where x = ^- (r + + r_ ) and y = 
p 2 + (z ± k) 2 and n = yjd + 5 
a = and 6 = we have, obviously, that the quadrupolar 
deformation (J7J) corresponds to the deformation for TS2 
solution, so all the consideration pointed out above are 
also valid for this case. 
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III. GEODESICS DYNAMICS FOR TEST 
PARTICLES 



Following an standard procedure, we define the La- 
grangian function C as 2C = g ^dx 11 / 'dr dx v / 'dr , for the 
stationary axisymmetric metric ((3]), we obtain 



k 2 



2C = f(t-^Y-j 

•2 



y 



i 



,2"; 



(x 2 -y 2 ) 



x 2 - 1 



2\ -2 

y w 



(11) 



where the overdot indicates derivation respect the proper 
time r. Using the Hamiltonian formalism we found that 
the motion equations for the test particle are given by 



fix 2 - 1) 



k 2 e 2 ^(a 



y 



-Px, y = - 



f(i-y 2 ) p 

k 2 e^{x 2 - y 2 ) v ' 
(12) 



P., 



E? 
f 
f 



f{L + Euf 



k 2 (x 2 - 1)(1 -y 2 ) 



P„ 



with 



fc 2 e 2 T(x 2 - y 2 ) 
1 (E 2 f{L + Eu) 2 

2\ f k 2 {x 2 ~l)(\-y 2 ) 
f 

k 2 e 2 ^(x 2 - y 2 ) 



[P 2 (* 2 -l)+P 2 (l-y 2 )] (13) 



[P 2 (x 2 -l)+P 2 (l-y 2 )]\(U) 



E = f(t-w<p) 
L = 



uf(i - ucp) - —(x 2 



(15) 

\){l~y 2 )^. (16) 



The constants of integration E and L are related to the 
energy and to the angular momentum of the test par- 
ticle, respectively. In the case of timelike geodesic, the 
lagrangian C satisfies C — 0.5, this relation allows us to 
define an effective potential, which explicitly is 



f 



k^e 2 -y{x 2 - y 2 ) 



E 2 
/ 



f(L + Em) 2 
fc 2 (x 2 -l)(l-y 2 ) 



1 



By the $ definition the motion must be restricted to the 
region $ > 0. With the aim of study the dynamic of 
geodesies is necessary to be sure that the test particles 
motion is in a confinement region. The existence of such 
regions is determined throughout the condition $ > 0. 



of the system, if we have E, L and x(0), 2/(0), P x (0), the 
momentum P y (0) will be determined for this equation. 
The values of x(0) and y(0) are selected in such form 
that satisfies the condition $ > for confined motions, 
then the only arbitrary parameter is P x (0). These con- 
stants of motion indicate to us that the geodesic motion 
is performed in a three dimensional effective phase space 
in which the Poincare section method is an adequate tool 
to study the motion (see Letelier & Gueron p7l[25|). 



Case I: TS2 solution 



First we shall to consider the TS2 solution. In this case, 
the k parameter introduced in f3| is completely deter- 
mined by the mass and the angular momentum of the 
source, implying that we only can vary the energy E, 
and the angular momentum L, of the test particle. For 
this solution, we only found bounded region of motion 
like that shown in Figflja). Any configuration with two 
or more bounded regions of motion was not found. In 
FigOJa), the curve $ = in the plane xy for E = 0.94, 
L = —3.12 is plotted and it is observed that there just ex- 
its one bounded region and one escape region. In FigOJb) 
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FIG. 1: (a) Boundary contour $ = using E = 0.94, L — 
—3.12. There is one escape zone in the left-hand side of the 
picture, which correspond to small values of x, and a closed 
zone of bounded motion to the right, (b) Poincare sections in 
the plane xP x for the values defined in (a). We have regular 
motion. 

we can see that the motion in the bounded region of the 
Figda) is completely regular such as that which occurs 
in the Schwarzschild, Kerr and Kerr-Newman black holes 
The geodesies for this solution were studied using 
surface sections for many different values of E and L. 
The numerical results suggest the existence of only inte- 
grable geodesies. 



The solution to the equation system p2 |) - ([T4|) . could be 
found using a symplectic Runge-Kutta method. That 
method let us find the numerical solution of the sys- 
tem, given the constants E, L and the initial conditions 
x(0), y(0), P x (0) and -P y (0). By the existence of the in- 
tegral of motion C = H = 0.5 with 7i the hamiltonian 



B. 



Case II: Manko et al. Solution 



In this case, we used the available information present in 
the literature for typical values for the multipolar struc- 
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ture of neutron stars. In particular, we took numerical 
data from the Berti and Stergioulas work [13]. In [13 ] . 
Berti and Stergioulas solved in a numerical way the full 
Einstein equations to determine the spacetime for rapidly 
rotating neutron star along of sequences of constant rest 
mass for selected equation of state (henceforth EOS) de- 
noted as EOS A Hzl, EOS AU [H, EOS FPS H, EOS 
L [H and EOS APRb 



(a) (b) 




FIG. 2: (a) Boundary contour $ = for M = 1.84OM , 
J = 3.683, b = -0.3792, E = 0.1 and L = -5.6. There is 
two escape zones, and a closed zone of bounded motion, (b) 
Poincare sections in the plane xP x for the values defined in 
(a). We see only regular motions. 
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FIG. 3: (a) Boundary contour $ = for M = 1.936M , 
J = 4.498, b = -0.3080, E = 0.95 and L = -9.0. There is 
one escape zone, and a closed zone of bounded motion, (b) 
Poincare sections in the plane xP x for the values defined in 
(a). The geodesies are only regular. 
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FIG. 4: (a) Boundary contour $ = for M = 1.936M©, 
J = 4.498, b = 0.8, E = 0.96 and L = -9.9. There are two 
small escape zones, and two closed zones of bounded motion, 
(b) Poincare sections in the plane xP x for the values defined 
in (a). We see chaotic motion to the left-hand side and regular 
motion to the right-hand side. 
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FIG. 5: (a) Boundary contour $ = for M = 1.936M©, 
J = 4.498, b = 0.8, E = 0.971 and L = -9.3. There are 
two small escape zones, and two regions linked by a narrow 
connection of bounded motion, (b) Poincare sections in the 
plane xP x for the values defined in (a). We see chaotic motion 
in the left-hand side of the figure and in a small external region 
on the right-hand side. 



Additionally, they matched the Manko et. al. solution 
@ to the numerical solutions imposing the condition that 
the mass-quadrupole moment of the numerical and ana- 
lytic spacetimes be the same. Under this condition, the b 
parameter is determined from the numerical data. From 
[l2l ] we can see that the Manko et al. solution fixes better 
for the EOS FPS, for that reason we took all the values 
presented by Berti and Stergioulas for each sequence of 
mass of the EOS FPS founding bounding potentials like 
the presented in FigH^a) and Figj3ja). 

In Figl2a), we present the boundary contour for $ = in 
the xy plane for M = 1.84OM , J = 3.683, b = -0.3792, 
E = 0.1 and L = —5.6, observing that exists two scape 
regions and a bounded region. In Fig(3]^a) boundary 
contour for $ = in the plane xy for M = 1.936A/ Q , 
J = 4.498, b = -0.3080, E = 0.95 and L = -9.0 is 
presented. This configuration is very similar to that pre- 
sented in the first case and is the commonest shape found. 
The geodesic behavior inside of the bounded region for 
these potentials is presented in the Figl^b) and Fig[3Jb), 
respectively. From there can be seen that the motion is 
completely regular. 

In order to found chaotic behavior, we let ourself mod- 
ify sightly the values of the b parameter changing the 
quadrupole deformation and also the properties of the 
test particle throughout the change in its energy and 
angular momentum values. Introducing this changes 
we found potential like the presented in FigfJJa) and 
FigEa). 

In FigHJa), we present the boundary contour for M = 
1.936M©, J = 4.498, b = 0.8, E = 0.96 and L = -9.9 ob- 
serving that there are two very small escape regions and 
two disconnected bounded regions. In Fig[5ja) bound- 
ary contour for $ = in the plane xy for M = 1.936A/©, 
J = 4.498, b = 0.8, E = 0.971 and L = -9.3 is presented. 
This configuration conserves the tiny scape regions but 
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the confined regions are connected. The geodesic behav- 
ior inside of the bounded region in the phase space for 
these potentials is presented in the FigQJb) and Fig[5^b), 
respectively. From there one can notice the presence of a 
mixed phase space -chaotic and regular- to the left but 
only regular to the right in both cases. 



IV. CONCLUDING REMARKS 



We have constructed the Poincare sections for a lot of 
possible combinations of energy and angular momentum 
(and deformation parameter in the Manko et al. case) 
which confine the motion of the test particles orbiting 
around two types of sources. 

We concluded that is apparently impossible to find 
chaotic geodesies around a source described for the TS 
solution and that the stability of the geodesies does not 
depend on any way of the relative spin direction of the 
center of attraction nor of the angular momentum of the 
test particle. This result shows numerical evidences of 
the existence of only integrable geodesies for this sys- 
tem. In other words, the case of the test particle turn- 
ing around a Tomimatsu-Sato source type is completely 
integrable. This work attempts to complete the study 
of the geodesic dynamics in the well known trilogy of 
axially symmetric solutions (Schwarzschild, Kerr, Tomi- 
matsu Sato 8 = 2). 

On the other hand, we use the Manko et al. solution 
with the idea of analyze the geodesic dynamics of test 
particles in rapidly rotating neutron stars with equation 
of state FPS. For this we took the values presented by 



Berti and Stergioulas for the parameters describing the 
multipolar structure of neutron stars. In this case we 
only found regular geodesies like in the precedent metric. 
In order to find chaotic behavior, we modify the values 
of the deformation parameter. Introducing these changes 
we found potentials with two bounded regions, in which 
the Poincare sections exhibit mixed phase spaces, i.e., 
phase spaces when some of the KAM tori survives inside 
a chaotic sea. 

The main result in the Letelier's work was the existence of 
chaotic geodesies in the geometry that characterizes the 
prolate case. In the oblate case these orbits appear to be 
regular. In our case we found chaotic geodesies (in the 
Manko's solution) only in the geometry that characterizes 
the oblate case. In [17| Letelier et al. analyzed the influ- 
ence of the introduction of multipolar terms correspond- 
ing to the quadrupole deformation in the Schwarzschild 
and Kerr solution finding chaotic behavior for some val- 
ues of the deformation. Considering the former state- 
ment and based on the numerical evidence of our work 
we could conjecture that the parameter of arbitrary de- 
formation is which what induces the ergodic motion for 
uncharged test particles orbiting general relativistic vac- 
uum sources. 
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